% This code evaluates firm j's FOC w.r.t. p_j1 at the candidate value x0.
function foc_out=foc_1st_lapse(x0, p2_jj, ...
    cost1_e_jj, cost2_1_jj, j, v_c, v_m, v0_mm, p1_target_jj,...
    D_claims2_jj, pr_s_row, ...
    alpha, Delta_xi_jj, ...
    KM, resource_y1, resource_y2, N_y, rho, crra, U0_mm, pr_y, sigma_c1, ...
    myopic, beta_c, beta_f, T1, B1_c, B2_c, B1_f, B2_f, delta, cost_weight, use_cost_weight)
%% firm j's new price
p1_jj = x0(1); % new initial price being solved for 
p2_jj_c = p2_jj; % for correct beliefs, fixed to previous iteration value
p2_jj_m = p1_jj*ones(1, KM); % for misbeliefs, updated to the new initial price

%% update demand using j's new initial price (unilateral deviation)
% consumption after paying the premium to j
c1 = resource_y1-p1_jj; % N_y x 1
c2_c = repmat(resource_y2,1,KM)-repmat(p2_jj_c,N_y,1); % N_y x KM 
c2_m = repmat(resource_y2,1,KM)-repmat(p2_jj_m,N_y,1); % N_y x KM 

force_lapse_c = (c2_c<1);
force_lapse_m = (c2_m<1);

% 2nd-period utility if keep contract woith j
U_stay_c = B2_c*(alpha*u_c(c2_c,rho,crra) +  Delta_xi_jj); % N_y x KM
U_stay_m = B2_c*(alpha*u_c(c2_m,rho,crra) +  Delta_xi_jj); % N_y x KM

% lapse utility: not dependent upon agg state
U_lapse = U0_mm;
U_lapse_mat = repmat(U_lapse, 1, 4); % N_y x KM

% expected utility in 2nd period
U_c = log(exp(U_stay_c) + exp(U_lapse_mat));% N_y x KM
U_c(force_lapse_c==1) = U_lapse_mat(force_lapse_c==1);

U_m = log(exp(U_stay_m) + exp(U_lapse_mat));% N_y x KM
U_m(force_lapse_m==1) = U_lapse_mat(force_lapse_m==1);

EU_c = U_c*pr_s_row'; % N_y x 1
EU_m = U_m*pr_s_row'; % N_y x 1

% choice-specific value for firm j
v_new_c = B1_c*(alpha*u_c(c1,rho,crra) +  Delta_xi_jj)+ (beta_c^T1)*EU_c; % N_y x 1
v_new_m = B1_c*(alpha*u_c(c1,rho,crra) +  Delta_xi_jj)+ (beta_c^T1)*EU_m; % N_y x 1

% fix consumers' values for other firms & update the value for firm j only
v_c(:,j) = v_new_c; % % N_y x NN, update the value for firm j
v_m(:,j) = v_new_m;

ind_temp = ones(1,size(v_c,2)); % 1 x NN
s1_ind_c = share_ind_ftn(v_c, ind_temp, v0_mm, sigma_c1); %N_y x (N_e+1)
s1_ind_m = share_ind_ftn(v_m, ind_temp, v0_mm, sigma_c1); %N_y x (N_e+1)

% N_y x (N_e+1)
s1_ind = (1-myopic)*s1_ind_c + myopic*s1_ind_m;
s1_ind_j_c  = s1_ind_c(:,j); % N_y x 1, demand for j when correct beliefs
s1_ind_j_m  = s1_ind_m(:,j); % N_y x 1, demand for j when incorrect beliefs
s1_ind_j    = s1_ind(:,j); % N_y x 1, integrated demand

s1_mkt_jj = sum(s1_ind_j.*repmat(pr_y, 1, 1))'; % 1 x 1. firm j's updated 1st period mkt share

%% evulate firm j's FOC w.r.t initial price
% for market share derivative, refer to supply_estimation/D_s_price.m
% for profit derivative, refer to supply_estimation/minloss_estimation.m
%---------------------
% D_s1_p1: 1 x 1.
% lapse+edit: misinformed consumers' wrongly perceived prob of keeping the
% contract
stay_m = exp(U_stay_m)./(exp(U_stay_m) + exp(U_lapse_mat)); % N_y x KM
agg_stay_m = sum(stay_m.*repmat(pr_s_row, N_y, 1),2);  % N_y x 1, integrateded over k

% D_v_p1: N_y x 1 (myopic consumers)
D_v_p1_c = -B1_c*alpha*D_u_c(c1, rho, crra)/sigma_c1;
% lapse+edit: plug in agg_stay_m 
D_v_p1_m = (-B1_c*alpha*D_u_c(c1, rho, crra)...
    - ((beta_c^T1)*B2_c*alpha)*D_u_c(c2_m(:,1), rho, crra).*agg_stay_m)/sigma_c1;

% D_s1_ind_p1: N_y x 1
D_s1_ind_p1_c = D_v_p1_c.*s1_ind_j_c.*(1-s1_ind_j_c);
D_s1_ind_p1_m = D_v_p1_m.*s1_ind_j_m.*(1-s1_ind_j_m);
D_s1_ind_p1 = myopic*D_s1_ind_p1_m + (1-myopic)*D_s1_ind_p1_c;

% D_s1_p1: 1 x 1.
D_s1_p1 = sum(D_s1_ind_p1.*pr_y)';

%---------------------
% D_s2_p1: 1 x KM. (k)th element = derivative of s_jk2 w.r.t. p_j1
% lapse+edit: lapse decisions always based on actual prices in t=2
stay_c = exp(U_stay_c)./(exp(U_stay_c) + exp(U_lapse_mat)); % N_y x KM
D_s2_p1 = sum(repmat(D_s1_ind_p1,1,KM).* stay_c.*repmat(pr_y, 1, KM)); % 1 x KM

% expected claim 
D_claim_av_jj=zeros(1, KM);
D_s2_ind_p1=repmat(D_s1_ind_p1,1,KM).* stay_c; % N_y x KM
for k=1:KM
    D_claim_av_jj(k) = sum(cost_weight.*repmat(D_claims2_jj(k), 1, N_y)'.*D_s2_ind_p1(:, k).*pr_y, 'all');
end

%---------------------
% profit derivative
D_profit1_p1 = B1_f*(s1_mkt_jj + p1_jj*D_s1_p1); % 1 x 1

delta_p = p2_jj-repmat(p1_jj,1,KM); % 1 x KM

if use_cost_weight==0
    D_profit2_p1_k = (p2_jj - D_claims2_jj).*D_s2_p1 + (delta_p>0).*cost2_1_jj.*delta_p; % 1 x KM
else
    D_profit2_p1_k = p2_jj.*D_s2_p1 - D_claim_av_jj + (delta_p>0).*cost2_1_jj.*delta_p; % 1 x KM
end

D_profit2_p1 = (beta_f^T1)*B2_f*sum(D_profit2_p1_k.*pr_s_row,2); % 1 x 1
D_profit = D_profit1_p1 + D_profit2_p1;
p1_deviation = p1_jj-p1_target_jj;

% FOC 
foc_out(1)= cost1_e_jj - (1/p1_deviation)*D_profit;





















